# include "scene.h"
# include "discr.h"
# include <stdio.h>
# include <stdlib.h>
# include <string.h>
# include <errno.h>
# include <math.h>

# ifdef _OMP_SOLVE_
# include <omp.h>
# endif


# ifdef SOLVE_ENERGY
extern double temp_advance(struct grid *g);
# endif

extern void init_cond(struct grid *g);
extern short print_result(unsigned int t_cnt, const struct grid *g);
extern unsigned int vel_advance(struct grid *g, unsigned int *pcnt_max/*, double *E_absr*/);
extern double inlet_velocity(double r, double t, short index);
extern double dens(double T);
extern double kvisc(double T);
extern double thm_cond(double T);

static short print_recov(const struct grid *g);
static FILE *file_open(FILE *lptr, const char *file, unsigned int t_cnt, unsigned int nfreq);

static double orlansky_outlet(struct grid *g, unsigned int t_cnt);
static double heat_balance(const struct grid *g, double E_absr);


extern short solve(struct grid *g)
 {
  unsigned int t_cnt, pcnt_max, pcnt_total;
  short sst_cnt, div_cnt;
  double pm_accum, pE_accum, E_absr;

# ifdef TRACE_FILE
  fprintf(g->trptr, "Entering into 'solve()'...\n");
# endif

  g->u = g->arr.q[0];
  g->v = g->arr.q[1];
  g->w = g->arr.q[2];
  g->p = g->arr.q[3];

# ifdef SOLVE_ENERGY
  g->T = g->arr.q[15];
  g->vhsrc_flag = VHSRC_FLAG;
# endif

  g->density0 = dens(T_INLET);
  g->kvisc0 = kvisc(T_INLET);
  g->Re0 = 0.5*W_MAX*d_NOZZLE/g->kvisc0;
  fprintf(g->lptr, "\nReynolds Number at Inlet, Re0 = %E\n", g->Re0); fflush(g->lptr);
  g->thmcond0 = thm_cond(T_INLET);

  init_cond(g);

/* =================== Calculation at different time-steps ==================== */

  pcnt_total = sst_cnt = div_cnt = 0;
  g->dptr = g->tptr = (FILE *)0;
  t_cnt = 0;

  while(g->t <= tMAX)
    {
     if(t_cnt%DFILE_FREQ == 0)
          {
           if(g->dptr) if(g->dptr != stderr) fclose(g->dptr);
           g->dptr = file_open(g->lptr, DISPLAY_FILE, t_cnt, DFILE_FREQ);
          }
     fprintf(g->dptr, "\n\n\t[[ t_cnt = %u, t = %E Sec ]]\n", t_cnt, g->t); fflush(g->dptr);

     if(t_cnt%IRFILE_FREQ == 0)
          {
           if(g->tptr) if(g->tptr != stderr) fclose(g->tptr);
           g->tptr = file_open(g->lptr, ittr_RES_FILE, t_cnt, IRFILE_FREQ);
          }
     fprintf(g->tptr, "\n\n[[ t_cnt = %u, t = %E Sec ]]; ", t_cnt, g->t); fflush(g->tptr);

     pm_accum = orlansky_outlet(g, t_cnt);
     fprintf(g->tptr, "pm_accum = %E%c; ", pm_accum, '%'); fflush(g->tptr);

     pcnt_total += vel_advance(g, &pcnt_max);

#    ifdef SOLVE_ENERGY
     E_absr = temp_advance(g);
     pE_accum = heat_balance((const struct grid *)g, E_absr);
     fprintf(g->tptr, "pE_accum = %E%c; ", pE_accum, '%'); fflush(g->tptr);
#    endif

     if(pcnt_max == POISS_CNT_MIN) ++sst_cnt; else sst_cnt = 0;
     if(pcnt_max == POISS_CNT_MAX+1) ++div_cnt; else div_cnt = 0;

     t_cnt++;
     g->t += g->dt;

  /* Printing result and recovery files. */
     if( t_cnt%RESF_FREQ == 0 || sst_cnt == SST_CNT_MAX ) if(print_result(t_cnt, g) < 0) exit(-1);

     if(pcnt_total >= RECOVF_POISS_CNT) if(pcnt_max <= POISS_CNT_MAX)
          {
           print_recov(g);
           pcnt_total = 0;
          }

#    ifdef SST_CNT_MAX
     if(sst_cnt == SST_CNT_MAX)
           {
            fprintf(g->tptr, "STEADY-STATE solution; "); fflush(g->tptr);
            fprintf(g->lptr, "\n\t[Solution reached steady-state]\n"); fflush(g->lptr);
            break;
           }
#    endif

#    ifdef DIV_CNT_MAX
     if(div_cnt == DIV_CNT_MAX)
           {
            fprintf(g->tptr, "DIVERGED solution; "); fflush(g->tptr);
            fprintf(g->lptr, "\n\t[Solution diverged\n"); fflush(g->lptr);
            break;
           }
#    endif
    }
/* --------------------------------------Time-loop ends------------------------------------------ */

  if(g->dptr) if(g->dptr != stderr) fclose(g->dptr);
  if(g->tptr) if(g->tptr != stderr) fclose(g->tptr);

# ifdef TRACE_FILE
  fprintf(g->trptr, "...\"solve()\" ends.\n\n");
# endif

  return g->t > tMAX ? 1:-1;
 }

//APPLY THE CONSTRAINS FROM Equation-2.5 OF THE ORLANSKY PAPAER IN THIS FUNCTION phase_vel() :-()
static double orlansky_outlet(struct grid *g, unsigned int t_cnt)
   {
    unsigned short i, j, k, ii, jj, kk, tn, nr_gap, nr, nth;
    double cdtdz, dt, *th, *r, q0, q1, m_inlet, m_outlet, m_accum;

#   ifdef TRACE_FILE
    fprintf(g->trptr, "Entering into 'orlansky_outlet()'...\n");
#   endif

    r = g->r; nr = g->nr; nr_gap = g->nr_gap; nth = g->nth; th = g->th;  dt = g->dt;

 /* Rearraning the memory to backup the physical quantities at current time-instant. */
    for(i = t_cnt%t_BACKUP, tn = k = 0; tn < t_BACKUP; tn++)
         {
          if(tn == i) k = t_BACKUP;
          g->outBC.u[tn] = g->memOut.q[  t_BACKUP+tn-i-k];
          g->outBC.v[tn] = g->memOut.q[2*t_BACKUP+tn-i-k];
          g->outBC.w[tn] = g->memOut.q[3*t_BACKUP+tn-i-k];
#         ifdef SOLVE_ENERGY
          g->outBC.T[tn] = g->memOut.q[4*t_BACKUP+tn-i-k];
#         endif
         }

 /* Backing up the quantities at current time-instant. */
    for(i = 0; i < nr_gap; i++) for(ii = i+NrG+nr-nr_gap, k = 0; k < zCells_BACKUP; k++) for(kk = k+NzG, jj = NthG, j = 0; j < nth; j++, jj++)
       {
        g->outBC.u[0][i][j][k] = g->u[ii][jj][kk-1];
        g->outBC.v[0][i][j][k] = g->v[ii][jj][kk-1];
        g->outBC.w[0][i][j][k] = g->w[ii][jj][kk-1];
#       ifdef SOLVE_ENERGY
        g->outBC.T[0][i][j][k] = g->T[ii][jj][kk-1];
#       endif
       }

#   ifndef RECOV_FILE
 /* At outlet, data at time "t-dt" and "t-2*dt" are kept same as at current time "t" when calculation starts without recovery file. */
    if(t_cnt == 0) for(tn = 1; tn < t_BACKUP; tn++)
    for(k = 0; k < zCells_BACKUP; k++) for(i = 0; i < nr_gap; i++) for(j = 0; j < nth; j++)
       {
        g->outBC.u[tn][i][j][k] = g->outBC.u[tn-1][i][j][k];
        g->outBC.v[tn][i][j][k] = g->outBC.v[tn-1][i][j][k];
        g->outBC.w[tn][i][j][k] = g->outBC.w[tn-1][i][j][k];
#       ifdef SOLVE_ENERGY
        g->outBC.T[tn][i][j][k] = g->outBC.T[tn-1][i][j][k];
#       endif
       }
#   endif

 /* Calculating quantities at the first adjacent ghost cells near the outlet at next time step using the Orlansky BC. */
    for(i = 0; i < nr_gap; i++) for(j = 0; j < nth; j++)
       {
     /* Calculating "tnext_ug[i][j]" at the outlet. */
    /*  cdtdz = (g->outBC.u[2][i][j][1] - g->outBC.u[0][i][j][1])/(g->outBC.u[0][i][j][1] + g->outBC.u[2][i][j][1] - 2*g->outBC.u[1][i][j][2]); */
        cdtdz = g->outBC.u[2][i][j][1] - g->outBC.u[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.u[0][i][j][1] + g->outBC.u[2][i][j][1] - 2*g->outBC.u[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_ug[i][j] = (1-cdtdz)*g->outBC.u[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.u[0][i][j][1]/(1+cdtdz);

     /* Calculating "tnext_vg[i][j]" at the outlet. */
     /* cdtdz = (g->outBC.v[2][i][j][1] - g->outBC.v[0][i][j][1])/(g->outBC.v[0][i][j][1] + g->outBC.v[2][i][j][1] - 2*g->outBC.v[1][i][j][2]); */
        cdtdz = g->outBC.v[2][i][j][1] - g->outBC.v[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.v[0][i][j][1] + g->outBC.v[2][i][j][1] - 2*g->outBC.v[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_vg[i][j] = (1-cdtdz)*g->outBC.v[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.v[0][i][j][1]/(1+cdtdz);

     /* Calculating "tnext_wg[i][j]" at the outlet. */
     /* cdtdz = (g->outBC.w[2][i][j][1] - g->outBC.w[0][i][j][1])/(g->outBC.w[0][i][j][1] + g->outBC.w[2][i][j][1] - 2*g->outBC.w[1][i][j][2]); */
        cdtdz = g->outBC.w[2][i][j][1] - g->outBC.w[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.w[0][i][j][1] + g->outBC.w[2][i][j][1] - 2*g->outBC.w[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_wg[i][j] = (1-cdtdz)*g->outBC.w[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.w[0][i][j][1]/(1+cdtdz);

     /* Calculating "tnext_Tg[i][j]" at the outlet. */
#       ifdef SOLVE_ENERGY
     /* cdtdz = (g->outBC.T[2][i][j][1] - g->outBC.T[0][i][j][1])/(g->outBC.T[0][i][j][1] + g->outBC.T[2][i][j][1] - 2*g->outBC.T[1][i][j][2]); */
        cdtdz = g->outBC.T[2][i][j][1] - g->outBC.T[0][i][j][1];
        if(cdtdz != 0) cdtdz /= g->outBC.T[0][i][j][1] + g->outBC.T[2][i][j][1] - 2*g->outBC.T[1][i][j][2];
        if(cdtdz > 1) cdtdz = 1;
        else if(cdtdz < 0) cdtdz = 0;
        g->outBC.tnext_Tg[i][j] = (1-cdtdz)*g->outBC.T[1][i][j][0]/(1+cdtdz) + 2*cdtdz*g->outBC.T[0][i][j][1]/(1+cdtdz);
#       endif
       }

 /* Mass coming into the domain through inlet between time t~t+dt. */
    m_inlet = 0.25*g->density0*r_NOZZLE*r_NOZZLE*W_MAX\
       *(dt-tFRAC*(exp(log(1-wFRAC)*(g->t+dt)/tFRAC) - exp(log(1-wFRAC)*g->t/tFRAC))/log(1-wFRAC))*(th[NthG+g->nth]-th[NthG]);

 /* Mass going out of the domain between time t~t+dt. */
    for(q0 = 0, i = 0, ii = i+NrG+nr-nr_gap; i < nr_gap; i++, ii++)
        {
         for(q1 = 0, j = 0, jj = NthG+j; j < nth; j++, jj++) q1 += (th[jj+1]-th[jj])*(g->outBC.w[0][i][j][0]+g->outBC.tnext_wg[i][j]);
         q0 += (r[ii+1]+r[ii])*(r[ii+1]-r[ii])*q1;
        }
    m_outlet = -0.25*g->density0*q0*dt;

    m_accum = m_inlet - m_outlet;

#   ifdef N_w_CORRECT
    if(t_cnt < N_w_CORRECT)
#   endif
       {
     // Correcting "w_outlet" so that "m_accum = 0" between time t~t+dt.
        q0 = 4*m_accum/(g->density0*dt*(r[NrG+nr]*r[NrG+nr]-r[NrG+nr-nr_gap]*r[NrG+nr-nr_gap])*(th[NthG+nth]-th[NthG]));
        for(i = 0; i < nr_gap; i++) for(j = 0; j < nth; j++) g->outBC.tnext_wg[i][j] -= q0;
       }

#   ifdef TRACE_FILE
    fprintf(g->trptr, "...'orlansky_outlet()' ends.\n");
#   endif

 /* Return % of dmdt-accumulation at next time-stage. */
    return m_accum*100/m_inlet;
   }



# ifdef SOLVE_ENERGY
static double heat_balance(const struct grid *g, double E_absr)
    {
     double E_inlet, E_outlet, E_vhsrc, E_heater, *th, *r;
     unsigned short i, j, ii, jj, nth;

#    ifdef TRACE_FILE
     fprintf(g->trptr, "Entering into 'heat_balance()'...\n");
#    endif

     nth = g->nth; th = g->th; r = g->r;

  /* Determining total energy created by the volume-heat-source. */
     E_vhsrc = (g->vhsrc_flag == 'y' ? \
                  (th[NthG+nth]-th[NthG])*d_CYL*d_CYL*h_CYL*g->density0*Cp_WATER*(EXPECTED_T_RISE/EXPECTED_TIME)*g->dt/8:0);

  /* Calculating heat input because of heater on the top surface. */
     E_heater = g->dt*TOP_HEATER_HF*0.5*g->r_heater*g->r_heater*(g->th[NthG+g->nth]-g->th[NthG]);

  /* Calculating the energy coming into the domain through the inlet in time t~t+dt. */
     E_inlet = 0.25*Cp_WATER*T_INLET*g->density0*r_NOZZLE*r_NOZZLE*W_MAX\
       *(g->dt-tFRAC*(exp(log(1-wFRAC)*(g->t+g->dt)/tFRAC) - exp(log(1-wFRAC)*g->t/tFRAC))/log(1-wFRAC))*(th[NthG+g->nth]-th[NthG]);

  /* Calculating the heat gone out of the domain through outlet between t~t+dt. */
     for(E_outlet = 0, i = 0, ii = NrG + g->nr - g->nr_gap; i < g->nr_gap; i++, ii++) for(j = 0, jj = NthG; j < g->nth; j++, jj++)
     E_outlet += -0.0625*(r[ii+1]+r[ii])*(r[ii+1]-r[ii])*(th[jj+1]-th[jj])*(g->outBC.w[0][i][j][0]+g->outBC.tnext_wg[i][j])\
                     *(g->outBC.T[0][i][j][0]+g->outBC.T[0][i][j][1]+g->outBC.tnext_Tg[i][j]+g->T[ii][jj][NzG]);
     E_outlet *= g->density0*Cp_WATER*g->dt;

#    ifdef TRACE_FILE
     fprintf(g->trptr, "...'heat_balance()' ends.\n");
#    endif

  /* Return % heat-accumulation in the domain between time t~t+dt. */
     return (1-(E_outlet+E_absr)/(E_inlet+E_heater+E_vhsrc))*100;
    }
# endif


static FILE *file_open(FILE *lptr, const char *file, unsigned int t_cnt, unsigned int nfreq)
    {
     FILE *ptr;
     char fname[MAX_FNAME_LEN];

     if(file == (char *)0) ptr = stderr;
     else{
          if(!strcmp(file, NULL_DEV)) strcpy(fname, NULL_DEV);
          else sprintf(fname, "%d.%s.%u_%u.dat", RUN_CNT, file, t_cnt, t_cnt+nfreq-1);

          if( (ptr = (FILE *)fopen(fname, "w")) == (FILE *)0)
                {
                 fprintf(lptr, "WARNING:: In 'file_open()': Unable to open '%s' for writing: %s.\nRedirecting output to 'stderr'.\n",\
                                        fname, strerror(errno));
                 ptr = stderr;
                }
         }

     return ptr;
    }


static short print_recov(const struct grid *g)
  {
   unsigned short i, j, k, errn, tn;
   FILE *ptr;
   char fname[MAX_FNAME_LEN];

#  ifdef TRACE_FILE
   fprintf(g->trptr, "Entering into 'print_recov()'...\n");
#  endif

   sprintf(fname, "%d.%s", RUN_CNT, W_RECOV_FILE);

   if((ptr = fopen(fname, "w")) == (FILE *)0)
     {
      errn = errno;

      fprintf(g->tptr, " {'%s' - recovery file writing error: %s}; ", fname, strerror(errn));
      fflush(g->tptr);

      fprintf(g->lptr, "WARNING:: In 'print_recov()': Unable to open the recovery-file '%s' for writing. No recovery file is being written: %s.\n", fname, strerror(errn));

      return -1;
     }


   fprintf(ptr, "# IMPORTANT: Don't edit this file. Below is: \"time-instant (in Sec)\"  \"time-step (last-->present) (in Sec)\"\n%0.15E    %0.15E\n", g->t, g->dt);
   fprintf(ptr, "\n# Variables: u(i,j,k, t), v(i,j,k, t), w(i,j,k, t), p(i,j,k, t)");
#  ifdef SOLVE_ENERGY
   fprintf(ptr, ", T(i,j,k, t)");
#  endif
   fprintf(ptr, "\n");

//    for(i = NrG; i < g->nr + NrG; i++) for(j = NthG; j < g->nth + NthG; j++) for(k = NzG; k < g->nz + NzG; k++)
   for(i = 0; i < g->nr + 2*NrG; i++) for(j = 0; j < g->nth + 2*NthG; j++) for(k = 0; k < g->nz + 2*NzG; k++)
       {
        fprintf(ptr, "%0.15E  %0.15E  %0.15E  %0.15E", g->u[i][j][k], g->v[i][j][k], g->w[i][j][k], g->p[i][j][k]);
#       ifdef SOLVE_ENERGY
        fprintf(ptr, "  %0.15E", g->T[i][j][k]);
#       endif
        fprintf(ptr, "\n");
       }

   for(tn = 0; tn <= t_BACKUP-2; tn++ )
       {
        fprintf(ptr, "\n# Data at outlet at time \"t - %d*dt\" (used to determine phase velocities while applying the Orlansky BC):: Variables: u(i,j,k, t), v(i,j,k, t), w(i,j,k, t)", tn+1);
#       ifdef SOLVE_ENERGY
        fprintf(ptr, ", T(i,j,k, t)");
#       endif
        fprintf(ptr, "\n");

        for(i = 0; i < g->nr_gap; i++) for(j = 0; j < g->nth; j++) for(k = 0; k < zCells_BACKUP; k++)
              {
               fprintf(ptr, "%0.15E  %0.15E  %0.15E", g->outBC.u[tn][i][j][k], g->outBC.v[tn][i][j][k], g->outBC.w[tn][i][j][k]);
#              ifdef SOLVE_ENERGY
               fprintf(ptr, "  %0.15E", g->outBC.T[tn][i][j][k]);
#              endif
               fprintf(ptr, "\n");
              }
       }

   fclose(ptr);

   fprintf(g->tptr, "(Recovery file '%s' - written); ", fname); fflush(g->tptr);

#  ifdef TRACE_FILE
   fprintf(g->trptr, "...'print_recov()' ends.\n");
#  endif

   return 1;
  }
